# -*- coding: utf-8 -*-
import matplotlib.pyplot as plt
from skimage import img_as_ubyte,img_as_float,io,color,exposure,util,restoration,measure
import numpy as np
from PIL import Image,ImageEnhance,ImageFilter
import matplotlib
from scipy import misc,ndimage


#对数变换----实现输入图像的颜色通道直方图（未进行对数变换之前）
#image_as_ubyte------图像像素转为0~255之间的整数
def plot_image(image,title=''):
    plt.title(title,size=8),plt.imshow(image)
    plt.axis('off')
def plot_hist(r,g,b,title=''):
    r,g,b=img_as_ubyte(r),img_as_ubyte(g),img_as_ubyte(b)
    plt.hist(np.array(r).ravel(),bins=256,range=(0,256),color='r',alpha=0.5)
    plt.hist(np.array(g).ravel(),bins=256,range=(0,256),color='g',alpha=0.5)
    plt.hist(np.array(b).ravel(),bins=256,range=(0,256),color='b',alpha=0.5)
    plt.xlabel('pixel value',size=20),plt.ylabel('frequency',size=20)
#
#im=Image.open('F:\\信本191莫巧丽\\图像处理\\2021-11-23\\images\\parrot.png')
##应用PIL模块中的point()
#im=im.point(lambda i:255*np.log(1+i/255))
##分离通道
#im_r,im_g,im_b=im.split()
#plt.style.use('ggplot')
#plt.figure(figsize=(15,5))
#plt.subplot(121),plot_image(im,'Original Image')
#plt.subplot(122),plot_hist(im_r,im_g,im_b,'histogram for RGB channels')
#plt.show()

#幂律变换
#im=Image.open('F:\\信本191莫巧丽\\图像处理\\2021-11-23\\images\\tiger.jpg')
#imm=img_as_float(np.array(im))
#设置y值
#gamma=5
#im1=imm**gamma
#plt.style.use('ggplot')
#plt.figure(figsize=(15,5))
#im[...,0]分离通道
#plt.subplot(121),plot_hist(imm[...,0],imm[...,1],imm[...,2],'histogram for RGB channels(input)')
#plt.subplot(122),plot_hist(im1[...,0],im1[...,1],im1[...,2],'histogram for RGB channels(output)')
#plt.subplot(121),plot_image(im,title='Origal')
#plt.subplot(122),plot_image(im1,title='after')A
#plt.show()

#对比拉伸
#使用PIL作为点操作
#加载rgb图像，划分不同颜色通道，并可视化不同颜色通道直方图
#im=Image.open('F:\\信本191莫巧丽\\图像处理\\2021-11-23\\images\\panda.png')
#使用PIL的point（）函数实现对比拉伸，变换函数由contrast()函数定义为分段线性函数
#def contrast(c):
#    return 0 if c<70 else(255 if c>150 else(255*c-22950)/48)
#
#im1=im.point(contrast)

#使用PIL中的ImageEnhance模块
#constras=ImageEnhance.Contrast(im)
#im1=np.reshape(np.array(constras.enhance(2).getdata()).astype(np.uint8),(im.height,im.width,3))
###划分通道A
##im_r,im_g,im_b=im1.split()
#im_r=im1[...,0]
#im_g=im1[...,1]
#im_b=im1[...,2]
#plt.style.use('ggplot')
#plt.figure(figsize=(15,5))
##调用上面定义的显示图片函数
#plt.subplot(121),plot_image(im1,'Image')
#plt.subplot(122),plot_hist(im_r,im_g,im_b,'histogram for RGB channels')
#
#plt.show()

#二值化
#1.固定阈值的二值化----PIL的point()函数
#im=Image.open('F:\\信本191莫巧丽\\图像处理\\2021-11-23\\images\\veg.jpg').convert('L')
#绘制像素直方图
#plt.hist(np.array(im).ravel(),bins=256,range=(0,256),color='g')
#plt.xlabel('Pixel values'),plt.ylabel('Frequency'),plt.title('Histogram of pixel values')
#plt.show()

#半色调二值化
#im=Image.fromarray(np.clip(im+np.random.randint(-128,128,(im.height,im.width)),0,255).astype(np.uint8))
#
#plt.figure(figsize=(12,12))
#plt.gray()
#plt.subplot(221),plot_image(im,'Aoriginal image'),plt.axis('off')
#th=[0,50,100,150,200]#设置阈值
#for i in range(2,5):
#    im1=im.point(lambda x:x>th[i],'1')#每个像素判断,输出模式，1-二值图
#    plt.subplot(2,2,i),plot_image(im1,'binary image eith threshold='+str(th[i]))
#plt.show()


#基于scikit-image的对比度拉伸和直方图均衡化
#灰度图
#img=color.rgb2gray(io.imread('F:\\python练习\\图像实战练习\\images\\whale.png'))
#使用进行直方图均衡化
#img_eq=exposure.equalize_hist(img)#全局
#img_adapteq=exposure.equalize_adapthist(img,clip_limit=0.03)#局部
##绘制
#plt.gray()
#images=[img,img_eq,img_adapteq]
#titles=['original','after histogram','after adapteq histogram']
#plt.figure(figsize=(15,10))
#for i in range(3):
#    plt.subplot(1,3,i+1),plot_image(images[i],titles[i])
#plt.figure(figsize=(13,5))
#for i in range(3):
#    plt.subplot(1,3,i+1),plt.hist(images[i].ravel(),color='g'),plt.title(titles[i],size=10)    
#plt.show()
    
#基于scikit-image的对比度拉伸和直方图均衡化比较
#matplotlib.rcParams['font.size']=8
#def plot_image_and_hist(image,axes,bins=256):
#    image=img_as_float(image)#转为float
#    axes_image,axes_hist=axes
#    axes_cdf=axes_hist.twinx()
#    axes_image.imshow(image,cmap=plt.cm.gray)
#    axes_image.set_axis_off()
#    axes_hist.hist(image.ravel(),bins=bins,histtype='step',color='black')
#    axes_hist.set_xlim(0,1)
#    axes_hist.set_xlabel('Pixel intentsity',size=10)
#    axes_hist.ticklabel_format(axis='y',style='scientific',scilimits=(0,0))
#    axes_hist.set_yticks([])
#    image_cdf,bins=exposure.cumulative_distribution(image,bins)
#    axes_cdf.plot(bins,image_cdf,'r')
#    axes_cdf.set_yticks([])
#    return axes_image,axes_hist,axes_cdf

#读取图像
#im=io.imread('F:\\python练习\\图像实战练习\\images\\pepper.jpg')
##对比拉伸
#im_rescale=exposure.rescale_intensity(im,in_range=(0,100),out_range=(0,255))
##全局
#im_eq=exposure.equalize_hist(im)
##局部
#im_adapteq=exposure.equalize_adapthist(im,clip_limit=0.03)
#绘制图像

#fig=plt.figure(figsize=(13,7))
#axes=np.zeros((2,4),dtype=np.object)
#axes[0,0]=fig.add_subplot(2,4,1)
#for i in range(1,4):
#    axes[0,i]=fig.add_subplot(2,4,i+1,sharex=axes[0,0],sharey=axes[0,0])
#for i in range(0,4):
#    axes[1,i]=fig.add_subplot(2,4,5+i)
#axes_image,axes_hist,axes_cdf=plot_image_and_hist(im,axes[:,0])
#axes_image.set_title('Low constrast image',size=14)
#y_min,y_max=axes_hist.get_ylim()
#axes_hist.set_ylabel('Number of pixels',size=14)
#axes_hist.set_yticks(np.linspace(0,y_max,5))
#axes_image,axes_hist,axes_cdf=plot_image_and_hist(im_rescale,axes[:,1])
#axes_image.set_title('Contrast imaage',size=14)
#
#axes_image,axes_hist,axes_cdf=plot_image_and_hist(im_eq,axes[:,2])
#axes_image.set_title('Histogram equalization imaage',size=14)
#
#axes_image,axes_hist,axes_cdf=plot_image_and_hist(im_rescale,axes[:,3])
#axes_image.set_title('adaptive equalization imaage',size=14)
#    
#axes_cdf.set_ylabel('Fraction of total intensity',size=14)
#axes_cdf.set_yticks(np.linspace(0,1,5))
#fig.tight_layout()
#plt.show()


#直方图匹配

def cdf(im):
    c,b=exposure.cumulative_distribution(im)
    c=np.insert(c,0,[0]*b[0])
    c=np.append(c,[1]*(255-b[-1]))
    return c

def hist_matching(c,c_t,im):
    pixels=np.arange(256)
    new_pixels=np.interp(c,c_t,pixels)
    im=(np.reshape(new_pixels[im.ravel()],im.shape)).astype(np.uint8)
    return im

#plt.gray()
#im=(color.rgb2gray(io.imread('F:\\python练习\\图像实战练习\\images\\whale.png'))*255).astype(np.uint8)#输入图像
#im_t=(color.rgb2gray(io.imread('F:\\python练习\\图像实战练习\\images\\lena.jpg'))*255).astype(np.uint8)#模板图像
#
#plt.figure(figsize=(11,9))
#plt.subplot(2,3,1),plot_image(im,'Input image')
#plt.subplot(2,3,2),plot_image(im_t,'Template image')
#c=cdf(im)
#c_t=cdf(im_t)
#plt.subplot(2,3,3)
#p=np.arange(256)
#plt.plot(p,c,'r.-',label='input')
#plt.plot(p,c_t,'b.-',label='template')
#plt.legend(prop={'size':15})
#plt.title('CDF',size=13)
#im=hist_matching(c,c_t,im)
#plt.subplot(2,3,4),plot_image(im,'out image with hist.Matching')#直方图匹配输出图像
#c1=cdf(im)
##绘制直方图匹配输出图像的累积分布函数
#plt.subplot(2,3,5)
#plt.plot(np.arange(256),c,'r.-',label='input')
#plt.plot(np.arange(256),c_t,'b.-',label='template')
#plt.plot(np.arange(256),c1,'b.-',label='output')
#plt.legend(prop={'size':14})
#plt.title('CDF',size=14)
#plt.show()    


#RGB图像的直方图匹配

#线性噪声平滑
#PIL平滑
#1.基于ImageFilter平滑
#i = 1
#plt.figure(figsize=(10,35))
#for prop_noise in np.linspace(0.05,0.3,6):
#    im = Image.open('F:\\python练习\\图像实战练习\\images\\mandrill.jpg')
#    # 在图像中选择5000个随机位置
#    n = int(im.width * im.height * prop_noise)
#    x, y = np.random.randint(0, im.width, n), np.random.randint(0, im.height, n)
#    for (x,y) in zip(x,y):
#        im.putpixel((x, y), ((0,0,0) if np.random.rand() < 0.5 else (255,255,255)))
#    #保存噪声图片
##    im.save('F:\\python练习\\图像实战练习\\images\\mandrill_spnoise_' + str(prop_noise) + '.jpg')
#    plt.subplot(6,2,i)
#    plt.imshow(im)
#    plt.title('Original Image with ' + str(int(100*prop_noise)) + '% added noise', size=15)
#    i += 1
    #BLUR模糊滤波
#    im1 = im.filter(ImageFilter.BLUR)
#    plt.subplot(6,2,i)
#    plt.imshow(im1)
#    plt.title('Blurred Image', size=20)
#    i += 1
#plt.show()

#2.基于盒模糊核均值化平滑
#加载上面保存的噪声为10%的图片
#im = Image.open('F:\\python练习\\图像实战练习\\images\\mandrill_spnoise_0.1.jpg')
#plt.figure(figsize=(14,7))
#plt.subplot(1,3,1)
#plt.imshow(im)
#plt.title('Original Image', size=30)
#plt.axis('off')
##n为核大小
#for n in [3,5]:
#    #设置盒模糊核
#    box_blur_kernel = np.reshape(np.ones(n*n),(n,n)) / (n*n)
#    #使用ImageFilter.Kernel()函数
#    im1 = im.filter(ImageFilter.Kernel((n,n), box_blur_kernel.flatten()))
#    plt.subplot(1,3,(2 if n==3 else 3))
#    plt.imshow(im1)
#    plt.title('Blurred with kernel size = ' + str(n) + 'x' + str(n), size=13)
#    plt.axis('off')
#plt.suptitle('PIL Mean Filter (Box Blur) with different Kernel size', size=13)
#plt.show()

#基于高斯模糊滤波器平滑
#加载上面保存的噪声为 20%的图片
#im = Image.open('F:\\python练习\\图像实战练习\\images\\mandrill_spnoise_0.2.jpg')
#plt.figure(figsize=(15,6))
#i=1
#for radius in range(1,4):
#    im1=im.filter(ImageFilter.GaussianBlur(radius))
#    plt.subplot(1,3,i),plot_image(im1,'radius='+str(round(radius,2)))
#    i+=1
#plt.suptitle('PIL Gaussian Blur with different Radius',size=15)
#plt.show()
    
#基于Scipy ndimage进行盒核与高斯核平滑比较
#im = Image.open('F:\\python练习\\图像实战练习\\images\\mandrill_spnoise_0.2.jpg')
#k=7#7*7的核
#im_box=ndimage.uniform_filter(im,size=(k,k,1))
#s=2
#t=(((k-1)/2)-0.5)/s
##sigma是高斯滤波核的标准差
#im_gaussian=ndimage.gaussian_filter(im,sigma=(s,s,0),truncate=t)
#fig=plt.figure(figsize=(14,8))
#plt.subplot(131),plot_image(im,'original image')
#plt.subplot(132),plot_image(im_box,'with the box filter')
#plt.subplot(133),plot_image(im_gaussian,'with the gaussian filter')
#plt.show()    
    
#非线性噪声平滑
#中值滤波器
#i = 1
#plt.figure(figsize=(13,14))
#for prop_noise in np.linspace(0.05,0.3,3):
#    im = Image.open('F:\\python练习\\图像实战练习\\images\\mandrill.jpg')
#    n = int(im.width * im.height * prop_noise)
#    x, y = np.random.randint(0, im.width, n), np.random.randint(0, im.height, n)
#    for (x,y) in zip(x,y):
#        im.putpixel((x, y), ((0,0,0) if np.random.rand() < 0.5 else (255,255,255)))
#    im.save('F:\\python练习\\图像实战练习\\images\\mandrill_spnoise_' + str(prop_noise) + '.jpg')
#    plt.subplot(6,4,i)
#    plot_image(im,'Original Image with '+ str(int(100*prop_noise)) + '% added noise')
#    i += 1
#    for sz in [3,7,11]:
#        im1 = im.filter(ImageFilter.MedianFilter(size=sz))
#        plt.subplot(6,4,i),plot_image(im1,'output(Median Filter size='+str(sz)+')')
#        i += 1    
#plt.show()
#    
#最大值滤波器和最小值滤波器
#im=Image.open('F:\\python练习\\图像实战练习\\images\\mandrill_spnoise_0.1.jpg')
#plt.subplot(1,3,1)
#plot_image(im,'Original Image with 10% added noise')
#sz=3#设置最大最小滤波器大小为3
#im1=im.filter(ImageFilter.MaxFilter(size=sz))
#plt.subplot(1,3,2),plot_image(im1,'Output(Max Filter size='+str(sz)+')')
#im1=im1.filter(ImageFilter.MinFilter(size=sz))
#plt.subplot(1,3,3),plot_image(im1,'Output(Min Filter size='+str(sz)+')')
#plt.show()
    
#scikit-image平滑（去噪）
#双边滤波器
#生成噪声图像
#im=color.rgb2gray(img_as_float(io.imread('F:\\python练习\\图像实战练习\\images\\hill.jpg')))
#sigma=0.155
#noise=util.random_noise(im,var=sigma**2) 
##plt.imshow(noise)  
##双边滤波去噪
#plt.figure(figsize=(14,12))
#i=1
#for sigma_sp in [5,10,20]:
#    for sigma_col in [0.1,0.25,5]:
#        plt.subplot(3,3,i)
#        plt.imshow(restoration.denoise_bilateral(noise,sigma_color=sigma_col,
#             sigma_spatial=sigma_sp,multichannel=False))
#        plt.title(r'$\sigma_r=$'+str(sigma_col)+r',$\sigma_s=$'+str(sigma_sp),size=16)
#        i+=1
#plt.show()
    
    
#非局部均值滤波器
#def plot_image_axes(image,axes,title):
#    axes.imshow(image)
#    axes.axis('off')
#    axes.set_title(title,size=20)
#parrot = img_as_float(io.imread('F:\\python练习\\图像实战练习\\images\\parrot.png'))
#
#sigma = 0.25
#noisy = parrot + sigma * np.random.standard_normal(parrot.shape)
#noisy = np.clip(noisy, 0, 1)
#
##估计噪声标准偏差从噪声图像
#sigma_est = np.mean(restoration.estimate_sigma(noisy, multichannel=True))
#
#patch_kw = dict(patch_size=5,      # 5x5 
#                patch_distance=6,  # 13x13 
#                multichannel=True)
# 非局部均值慢
#denoise = restoration.denoise_nl_means(noisy, h=1.15 * sigma_est, fast_mode=False, **patch_kw)
## 非局部均值块
#denoise_fast = restoration.denoise_nl_means(noisy, h=0.8 * sigma_est, fast_mode=True, **patch_kw)
#
#fig, axes = plt.subplots(nrows=2, ncols=2, figsize=(10, 7), sharex=True, sharey=True)
#plot_image_axes(noisy,axes[0,0],'noisy')
#plot_image_axes(denoise,axes[0,1],'non-local means\n(slow)')
#plot_image_axes(parrot,axes[1,0],'orignal\n(noise free)')
#plot_image_axes(denoise_fast,axes[1,1],'non-local means\n(fast)')
#
#fig.tight_layout()
#psnr_noisy =measure.compare_psnr(parrot, noisy)
#psnr = measure.compare_psnr(parrot, denoise.astype(np.float64))
#psnr_fast = measure.compare_psnr(parrot, denoise_fast.astype(np.float64))

#信噪比
#print("PSNR (noisy) = {:0.2f}".format(psnr_noisy))
#print("PSNR (slow) = {:0.2f}".format(psnr))
#print("PSNR (fast) = {:0.2f}".format(psnr_fast))

#plt.show()    
    
#Scipy ndimage 平滑
#lena=io.imread('F:\\python练习\\图像实战练习\\images\\lena.jpg')
##添加随机噪声
#noise=np.random.random(lena.shape)
#lena[noise>0.9]=255
#lena[noise<0.1]=0
#plot_image(lena,'noise image')
##plt.show()
#fig=plt.figure(figsize=(11,8))
#i=1
#for p in range(25,100,25):
#    for k in range(5,25,5):
#        plt.subplot(3,4,i)
#        #使用percentile_filter平滑
#        filtered=ndimage.percentile_filter(lena,percentile=p,size=(k,k,1))
#        plot_image(filtered,str(p)+'percentile,'+str(k)+'x'+str(k)+'kernel')
#        i+=1
#plt.show()
#    
    
    
    
    